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ABSTRACT 

A near-equal mass binary black hole can clear a central cavity in a circumbinary 
accretion disk; however, previous works have revealed accretion streams entering this 
cavity. Here we use 2D hydrodynamical simulations to study the accretion streams 
and their periodic behavior. In particular, we perform a suite of simulations, covering 
different binary mass ratios q — M2/M1 in the range 0.01 ^ q ^ 1. In each case, 
we follow the system for several thousand binary orbits, until it relaxes to a stable 
accretion pattern. We find the following results: (i) while the binary is efficient in 
maintaining a low-density cavity, the time-averaged mass accretion rate into the cavity, 
through narrow coherent accretion streams, is suppressed by at most a factor of ~ 5, 
compared to a disk with a single BH with the same mass; (ii) the largest suppression 
occurs for q « 0.05; binaries whose mass ratios are either lower or higher both suppress 
accretion less significantly; (iii) for q > 0.05, the accretion rate is strongly modulated 
by the binary, and depending on the precise value of q, the power spectrum of the 
accretion rate shows either one, two, or three distinct periods; and (v) for q < 0.05, 
the accretion rate becomes steady, with no time- variations. Most binaries produced 
in galactic mergers are expected to have q > 0.05. If the luminosity of these binaries 
tracks their accretion rate, then a pcriodogram of their light-curve could help in their 
identification, and to constrain their mass ratio and disk properties. 



Key words: black hole physics 
gravitational waves 



accretion, accretion discs — galaxies: active 



1 INTRODUCTION 

Massive black holes (MBHs) appear to reside in the nuclei 
of most nearby galaxies (see, e.g., reviews by Kormendy & 
|Richstone|1995| and |Ferrarese fe~F ord 2005). In hierarchical 
structure formation models, galaxies are built up by merg- 
ers between lower-mass progenitors, which deliver nuclear 
MBHs (e.g. |Springel et ~aL] |2005"1 |Robertson eTaL] [2006]) , 
along with a significant amount of gas ( Barnes & Hernquist 



1992), to the central region of the newly born post-merger 



resolvable orbital separations. Another possible hindrance, 
which we address in this paper, is that the outward gravita- 
tional torques from the binary can balance the inward vis- 
cous and pressure torques, clearing a central cavity in a pu- 
tative circumbinary gas disk ( Artymowicz fc Lubow|[l994 1, 
possibly rendering the system too dim for detection. Overall, 
identifying MBHBs is difficult, and a better understanding 
of their expected observational signatures, especially those 
based on time- variability (Haiman et al. 20091, is needed. 



galaxy. Since mergers are common (e.g. |Haehnelt fc Kaufi- 
mann 2002 ), it follows that massive black hole binaries (MB- 
HBs) should also be common in galactic nuclei. 

Despite this expectation, observational evidence for 
MBHBs remains scarce (see, e.g., a review by Komossa 



Merging MBHBs should be unambiguously identifiable by 



gravitational wave (GW) detectors, such as eLISA ( Amaro- 
Seoane fc et. al|2012 1 or ongoing Pulsar Timing Arrays (e.g. 



2006). The dearth of MBHBs could be attributed to sev- 



eral factors: it is possible that typically only one of the two 
BHs is active at spatially resolvable separations; binaries 
may also lose their angular momentum efficiently due to the 
surrounding stars and gas and quickly move to spatially un- 
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Lommen|2012[ ). Identifying the electromagnetic (EM) coun- 
terparts of these GW sources (among the many false candi- 
date galaxies in the GW error box) will, however, likewise 
require an understanding of their observational signatures. 

Recent studies have explored the gasdynamics of cir- 
cumbinary accretion disks around near-equal mass binaries 
in some detail. Since the system is not axisymmetric, this 
requires a two- or three-dimensional treatment. Mac Fadyen| 
|fc Milosavljevic| ( |2008| ) (hereafter MM08) have run two- 
dimensional hydrodynamical simulations for an equal mass 
binary. Three-dimensional smoothed particle hydrodynam- 
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ical (SPH) simulations have been carried out for equal-mass 
and 2:1 mass-ratio binaries by Hayasaki et al. (|2007| ) and 
for a 3:1 mass-ratio binary by Cuadra et al. ( |2009| ) and 
Roedig et aL| (|2012[). |Shi et al'| (|2011[) have followed up 



on the work of MM08 for an equal-mass binary by run- 
ning 3D magneto-hydrodynamical (MHD) simulations, and 
Noble et al. ( 2012[ ) have further added a post- Newtonian 
treatment of general relativistic (GR) effects, and followed 
the disk through the late stages of orbital inspiral (from 
orbital separation r = 20M to 8M).|Farris et al.|([20TT| fol- 



lowed the merger of an equal-mass binary and a surround- 
ing disk through merger in full 3D general relativistic MHD, 



starting from 10M. Finally, Farris et al. (20121 have added 



gas-cooling to GRMHD simulations of an equal-mass binary 
prior to decoupling, through decoupling and to merger start- 
ing from 10M. A generic result of all of these studies is that 
a low-density cavity is carved out by the binary torque, but 
gas leaks into the cavity through non-axisymmetric streams. 
These streams can power significant accretion onto the bi- 
nary components, and should lead to bright EM emission. 

A particularly promising feature is that the accretion 
rate onto the BHs can be both high and strongly variable, 
modulated by the binary's orbital motion. This could allow 
a detection of sub-pc binaries by looking for periodic varia- 
tions in the luminosity of AGN-like objects ( jHaiman et al. 
2009) or periodic shifts and intensity variations of spectral 
lines (e.g. Haiman et al.|2 009 She n fc Loeb| 2010 Eracleous 



et al.|[2011| and references therein). If the accretion remains 
significant and periodic down to <Cpc separations, then it 
could also enable the identification of EM counterparts of 
gravitational wave sources: either for p recursors to eLISA 
1O 7 M range ( |Kocsis et al.|2006 



sources in the M — 10 J 



2008) or by detecting periodic modulations of more mas- 
1O 9 M binaries discovered by pulsar timing 



sive M = 10' 

arrays (PTAs; [Taiiaka et al.|2012| |Sesana et al.||2012 1 

Although existing studies have focused on near-equal 
mass MBHBs, in reality, coalescing MBHBs should have 
a distribution of mass ratios q = M2/M1Q Mergers oc- 
cur between galaxies over a wide range of sizes, harboring 
central BHs of different masses, so that MBHBs resulting 
from galactic mergers should have a correspondingly wide 
range of mass ratios. Studies based on Monte-Carlo realiza- 
tions of dark matter merger trees indeed find broad distri- 
butions between 10 -2 < q < 1, genera lly peaking in the 
Volonteri et al. 2003 Sesana et al 



range q ~ 0.1 — 1 (e.g. 



2005 2012 Gergely fc Biermann 20121. However, the pre- 



dictions depend on the occupation fraction of MBHs, the 
redshift-evolution of the correlation between the masses of 
MBHs and their host galaxies, as well as on the limit on the 
mass ratio of host galaxies whose nuclear MBHs can coa- 
lesce; q < 0.1 mergers could in fact be most common (e.g. 
Lippai et al.|2009| |. 

Here we follow up on the earlier work of MM08 and 
move beyond the near-equal mass binary case. We study the 
periodicity and the time-averaged rate of accretion across 
the central cavity, by running 2D hydrodynamical simula- 
tions of a circumbinary disk for 12 different binary mass 
ratios ranging from q — 0.01 to q = 1. Clearly, one expects 



that in the limit q — > 0, the accretion rate approaches that of 
an accretion disk around a single BH, and will no longer be 
time-variable. The main goal in this paper is to answer the 
following basic questions: How does the mean accretion rate, 
and its fluctuations, depend on the mass ratio ? In particu- 
lar, down to what mass ratio is the mean accretion and/or 
its variability significantly affected by the binary torques? 

The rest of this paper is organized as follows: In § [2] we 
describe the setup of our numerical simulations, including a 
few changes we made to the public version of the Eulerian 
grid code Flash, basic tests of the modified code, and the 
initial and boundary conditions we adopted. In § [3j we dis- 
cuss scaling the simulations to physical parameters, such as 
black hole mass and orbital separation, and discuss the cor- 
responding orbital and residence times. In § [4] we present 
our main results, where we find four distinct patterns for 
the time-variability of the accretion rate, as a function of 
the mass ratio q. We also discuss the interpretation of these 
results, as well as some caveats, in this section. Finally, in 
§[5] we conclude by briefly summarizing our main results and 
their implications. 



2 DETAILS OF NUMERICAL SIMULATIONS 

To simulate a gas disk in the gravitational field of a bi- 
nary, we use the Eu lerian grid-based h ydrodynamical code 
FLASH (Version 3.2; |Fryxell et al.|2000[ )PlFLASH solves the 
volume-integrated fluid equations by solving the Riemann 
'shock-tube' problem at each cell boundary. A piece-wise 
parabolic representation of the fluid variables is used to in- 
terpolate between cells, i.e. FLASH is a PPM code, accurate 
to 2nd order in both space and time. FLASH uses a mono- 
tonicity constraint, rather than artificial viscosity, to control 
oscillations near discontinuities. This makes it well suited for 
following supersonic fluid dynamics in the inner regions of 
circumbinary disk. FLASH also supports polar coordinates, 
which is convenient for simulating disks. 



2.1 Numerical Implementation and Assumptions 

We assume a geometrically thin accretion flow onto the bi- 
nary. This permits a decoupling of the fluid equations in the 
z direction, perpendicular to the plane of the disk, so that we 
can define height-integrated fluid variables and set up sim- 
ulations in two dimensions. In follow-up studies, we plan to 
extend these simulations to the full three dimensions, which 
we expect will be important in determining the amount of 
inflow into a putative circum-binary cavity (for recent 3D 
grid-based simulations, see |Shi et al.|[20TT| and [Noble et al.| 
2012|. In the present study, we choose 2D polar coordinates 
(r, <j>) and employ FLASH to solve the following standard 



1 |Nixon et al.| | |201l| explored a range of mass ratios, using 3D 
SPH simulations, but restricted their study to retrograde disks. 



2 We note that MM08 used an earlier release, Version 2, of the 
same code. 
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set of 2D hydrodynamical equations: 
<9E 

at 

dv I 



+ V ■ (Eu) = 



at 

d(ZE) 
dt 



+ (v- V)« 



VP - V$bin + V ■ vSIv 



+ V ■ [(E£ + P)v\ = Zv- (-V$ bi , 



(1) 



Here E is the vertically integrated disk surface density, 
v = Vrf + v^cj) is the fluid velocity, P is the pressure, v is the 
coefficient of kinematic viscosity, E is the total internal plus 
kinetic energy of the fluid, E = e + §ji?j 2 , and $bin is the 
gravitational potential of the binary. The gravitational po- 
tential is inserted into the simulation by hand, and is given 
by 



$btn( 



r, 0) = — ■ 



GM(1 + q)~ 



r 2 + 



-,1/2 



l+q- 1 J 1+q 

GMil + q- 1 )- 1 



2ra TCOs(0-fW) 



1+9 



+ COS (0 - fibint) 



1/2' 



(2) 



Here n bin = (GM/a 3 ) 1/2 is the binary's orbital frequency, 
a is the separation of the binary, M v > M s are the masses of 
the primary and the secondary, M — M p + M s is the total 
mass, and q = M s /M p ^ 1 is the mass ratio. The origin of 
the coordinate system is chosen coincide with the binary's 
center of mass. Note that we do not evolve the orbital pa- 
rameters of the binary (this assumption is justified for our 
physical parameter choices; see discussion in §[3] below). 

We neglect self-gravity of the disk. Given the local 
sound speed c s , the Toomre parameter Q = c s Q/ (nCS) 
can be written as Q ~ (H/r)(M/Md), for a disk with mass 
Md and vertical scale-height H in hydrostatic equilibrium. 
For a thin disk H/r < 0.1, but in all of our simulations, we 
choose M 2> Md and thus Q > 1, making the disk stable 
to gravitational fragmentation. This is justified for standard 
Shakura-Sunyaev disks, when the simulations are scaled to 
physical BH masses and sufficiently small separations (see 
§[3] bel ow). For comparison, we note that in their SPH simu- 
lations, [Cuadra et al.| ( [2009| ) and |Roedig et al.| ( |2012| ) studied 
more massive disks, with Md ~ 0.4M, making self-gravity 
important. 

The pressure is given, as in MM08, by a locally isother- 
mal equation of state, 



P = c 2 (r)E 



(3) 



where the sound speed is assumed to scale with radius as 
c a oc r- -1//2 . For a Keplerian potential, this corresponds to a 
constant disk scale-height to radius ratio, H/r = c s /v$ = 
M^ 1 , where v^, is the orbital velocity in the disk and M 
is the corresponding Mach number. Throughout our sim- 
ulations, we chose the disk sound speed such that, for a 
Keplerian azimuthal velocity, H/r = 0.1 (or M = 10) ev- 
erywhere]^] 



3 Since we actually add the binary's quadrupole potential, and 
the non-zero pressure makes the azimuthal velocities slightly sub- 
Keplerian, in practice Ai in the simulations approaches S3 18 near 
r = a and becomes a constant Ai ~ 10 at r > 5a. 



To incorporate viscosity, FLASH calculates the momen- 
tum flux across cell boundaries due to viscous dissipation 
(the last term in the second of equations [T|. To compute v 
we adopt the a prescription, v = ac s H (Shakura & Sunyaev 



19731, where a is a dimensionless parameter indicating the 
scale of turbulent cells, and the scale height is computed 
from H = c s r/v^ = O.lr with v$ being the Keplerian value. 
Following MM08, we chose a constant a = 0.01 (but see 
discussion in § 2.4 below). The viscosity implementation in 
FLASH3 follows the method of |Edgar| pOTJri] ), and is very 
different from that of FLASH2. Since neither implementa- 
tion is fully supported in polar coordinates, we made adjust- 
ments to the routines that compute the momentum flux and 
viscous diffusion time from u, originally in Cartesian coor- 
dinates. These changes amount to keeping track of volume 
elements when taking azimuthal derivatives. The method is 
outlined in polar coordinates in Edgar (20061. We test our 
viscosity implementation in § |2.4| below. 



2.2 Numerical Parameter Choices 

Unless specified otherwise, the inner edge of the computa- 
tional domain is chosen at r m i n = a and the outer edge at 
r max = 100a. Although we are only interested in the inner 
few r/a of the disk, extending the computations to larger 
radii acts as a buffer for small initial numerical transients, 
and also provides the reservoir of gas from which the inner 
regions can be fed. As in other Eulerian codes, FLASH uses 
a boundary zone of 'guard cells' which enforces boundary 
conditions. We chose a diode-type inner boundary condition: 
the values of fluid variables in cells bordering the guard cells 
are copied into the guard cells, with the restriction that no 
fluid enter the domain, v r (r m in) ^ 0. We adopt an 'outflow', 
outer boundary condition: this is identical to the diode ver- 
sion, except that flow is allowed both into and out of the 
domain. In practice, with our initial conditions in § |2.3| we 
find an outflow at the outer boundary of our disk. However, 
the simulation does not run for a significant fraction of a 
viscous time at the outer radius (see below) and we expect 
inflow at the outer boundary to establish itself eventually. 

The spatial resolution was fixed throughout the 
grid, with the default values set at [Ar/a, A<f)/2n] ~ 
[0.024, 0.0078], corresponding to a grid of » 4096 x 128 cells. 
The wavelength of a density wave due to a Lindblad reso- 
nance is of order 2ty H ~ 0.6r (e.g. |Dong et al.|2011| |Duffell 
fc MacFadyen|2012| >; as in MM08, the radial resolution was 
chosen to resolve this length scale with many cells. To test 
numerical convergence (see 



4.3.3 below), we performed sev- 



eral additional runs, increasing the radial and azimuthal res- 
olutions, by factors of 1.42 and 1.42 2 w 2.0, from the lowest 
resolution runs. The time resolution is set to be half of the 
shortest propagation time (viscous or dynamical) across a 
cell. The effects of changing the resolution will be discussed 
in § [4X3] below. 

We run our simulations for between 4 x 10 3 and 10 4 
binary orbits. For reference, we note that the viscous time 
can be related to the orbital time as 



_ 2r 2 _ M 2 

3 V 37TQ 



1060 



100 ) 



0.01 

a 



(4) 



Thus, our typical run of 4000 orbits corresponds to 
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Table 1. Summary of our simulation runs. Low, medium, and high radial and azimuthal resolutions are denoted by "LoAr"; "MidAr"; 
"HiAr" and "LoAtf>"; "MidA<j>"; "HiA<f>" and correspond to Ar/a = 0.035,0.24,0.017 and A^/2vr = 0.0078,0.0052,0.0039, respectively. 
The viscosity parameter a is set to 0.01 unless otherwise specified. 



Mass ratio q 


Spatial Resolution [Ar, A(p] 


No. Orbits (N olb = i max fi bin /27r) 


1.0 


[LoAr, LoA0], [MidAr, LoA</>] 


5000, 5000 


1.0 (a = 0.05) 


[LoAr, LoA0], [MidAr, LoA0] 


4500, 4500 


0.75 


[LoAr, LoA0], [MidAr, LoA0] 


4500, 4500 


0.5 


[LoAr, LoA0], [MidAr, LoA^] 


4500, 4500 


0.25 


[LoAr, LoA</>], [MidAr, LoA<£] 


4800 


0.1 


[LoAr, LoA0], [MidAr, LoA^j, [MidAr, Mid A<f>], [HiAr, Hi A<f>], 


4500 , 4500, 4500, 4200 


0.075 


[LoAr, LoA0], [MidAr, LoA^j 


4500, 4500 


0.065 


[LoAr, LoA0] 


4500 


0.06 


[LoAr, LoA0] 


8000 


0.05 


[LoAr, LoA</>], [MidAr, LoA<£] 


4500, 5000 


0.04 


[MidAr, LoA</>] 


4500 


0.03 


[MidAr, LoA</>] 


4500 


0.01 


[LoAr, LoA0], [MidAr, LoA<£] 


7000, 5000 





[LoAr,LoA</>], [MidAr, LoA</>], [MidAr, Mid Atj>], [HiAr, HiA</>], 


7000 , 5000, 4500, 4200 



viscous times at the innermost regions, but less than one 
viscous time at r > 3a. 

We have performed 28 runs altogether, for 12 different 
binary mass ratios between q = 0.01 and g = 1.0 (as well 
as control runs for g = 0, i.e. a single BH). The mass ratio, 
resolution, and the number of binary orbits followed in each 
of our simulation runs are summarized in Table [T] 



2.3 Initial Conditions 

In our initial conditions, we insert a central cavity around 
a binary, and also include a density pile-up just outside the 
cavity wall. Such a surface density profile is expected to 
develop during the inward migration of the secondary. When 
the secondary arrives at the radius where the local disk mass 
is too small to absorb the secondary's angular momentum, 
its migration stalls, the inner disk drains onto the primary, 
and continued accretion from larger radii causes a pile-up 



of gas outside the secondary's orbit (|Syer & Clarke| 1995| 


Ivanov et al.||1999 


|Milosavljevic & Phinney||2005| 


Chang 


et al. 


2010 Rafikov 


2012 1. As emphasized recently by 


Kocsis 


et al. 


||2012a|), the details of this process are still uncertain, 



as the coupled time-dependent migration, cavity formation, 
and pile-up, has not been modeled self-consistently, even in 
one-dimensional calculations. However, using self-consistent 



steady-state solutions, Kocsis et al. ( 2012b I showed that in 



many cases, the pile-up can cause overflow already at large 
binary separations. 

For simplicity, and for ease of comparison, we adopt 
the same initial surface density profile as in MM08. This 
profile was motivated by the earlier results of Milosa vljevic| 
& Phinney ( 2005 I; it has very little gas inside of r ~ 3a, and 



peaks at ~ 8a, 



S(r,t ) = So (y) 3 exp 



Here So is an arbitrary constant, and r s = 10a. The ini- 
tial profile is shown by the [blue] dashed curves in Figure [4] 
below. 

We also follow MM08, and incorporate pressure gradi- 
ents and the quadrupole contribution of the binary's poten- 



tial into the initial azimuthal velocity, 
3 



1 + 



(l + <?) 2 



1 dP (K\ 

+ ^dr~ (6) 



and account for viscous drift in the initial radial velocities, 



v r 



d / 3 dfi\ 
dr I dr J 



rS- (r 2 Q) 
dr v ' 



(7) 



We emphasize that with this profile, the disk is not 
initially in equilibrium, and material diffuses away from the 
peak of the surface density, both inward and outward, due to 
pressure gradients. However, after running the simulations 
for several thousand orbits, the system relaxes to a steady 
pattern of accretion, and we do not expect the initial profile 
to significantly influence our conclusions. 



2.4 Tests of Code Implementation 

To ensure that non-axisymmetric perturbations are not 
induced artificially in the disk, we simulate a standard 
Shakura-Sunyaev a-disk around a single black hole for 
10 4 binary orbits. To keep axisymmetry during these 
runs, we re-implemented the FLASH2 routine UNBIASED- 
GEOMETRY (UBCQ into the FLASH3 routines. UBG 
cleans up round-off errors in the cell boundary positions, 
and forces the polar grid to keep cell sizes uniform in the 
azimuthal direction. Without UBG, small perturbations in 
the azimuthal direction grow to significant size after a few 
hundred orbits. We have verified that after re-implementing 
UBG, axisymmetry is preserved for 10 4 binary orbits. 

To test our implementation of viscosity in polar coor- 
dinates, we ran simulations of a viscously spreading ring of 
matter in a Keplerian potential. When the kinematic vis- 
cosity coefficient v is constant in space and time, there ex- 
ist analytic solutions for the evolution of the surface den- 



sity S(a;,T) and radial velocity v r (x, T) ( |Pringle|198Tj ). The 
dimensionless radial and temporal variables x = r/ro and 
T = Ylvtr^ 2 , are scaled to the initial position of the ring 
ro, and to the viscous time at vq, respectively. 



4 This routine is titled cleanAastJiits in the FLASH2 download. 
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Figure 1. Tests of our viscosity implementation in FLASH3.2, 
for a viscously spreading ring of matter in a Kcplcrian potential. 
The solid curves represent analytic solutions ([Pringlc 198111 and 
the dashed curves show the results of numerical simulations, at 
four different times, as labeled. T = Ylvtr^ is the dimensionless 
time, scaled to the viscous time at the initial position ro of the 
ring. The top panel shows the surface density (in arbitrary units) 
as a function of the dimensionless radius x = r/ro. The bottom 
panel shows the velocity profile in units of x/T. 



To test our code, we replaced the potential of the bi- 
nary (equation [5| by that of a point mass, and initialized 
FLASH with Keplerian azimuthal velocity, surface density 
~E(x,Tq = 0.32), radial velocity v r (x,To = 0.32), a constant 
temperature such that v^/c a ~ 30 so the disk is thin, and 
a large, constant kinematic coefficient of viscosity such that 
the viscous time is of order an orbital time at the disk inner 
edge t v i ac r2if (V m in)/27r — 1. The latter choice of a large vis- 
cosity coefficient amplifies the effects of viscosity, and also 
helps avoid a small initial transient which emanates from 
the inner edge of the disk after ~ 1 orbital period. If this 
transient perturbs the initial conditions, agreement with the 
analytic solutions is more difficult to gauge. 

Figure [I] displays the results of the test. Our implemen- 
tation is in good agreement with the analytic solutions for 
the entire viscous timescale shown in the figure. However, 
we did find a peculiarity: this good agreement is achieved 
by choosing a value of v in the simulations that is a fac- 
tor of ~ 3 greater than used in the analytic solution. We 
suspect that this factor of 3 discrepancy is caused by the 
different FLASH2 and FLASH3 viscosity implementations, 
but we were unable to track down its origin. We also note 



that the similar "fudge factor" of ~ 2 — 3 is needed for our 
equal-mass binary runs to agree with the accretion rates and 
surface densities observed in MM08. We therefore trust the 
results of the above tests, and simply adopt coefficients of 
kinematic viscosity which are 3 times larger than the desired 
input value (i.e., to simulate a disk with a = 0.01, we in- 
put a = 0.03 into FLASH3). In section 



4.3.2 



we further 

discuss the effects on our results of changing the viscosity 
coefficient. 

Finally, we note that the above test does not exercise 
any azimuthal components of the viscous momentum flux, 
which are present in a non-axisymmetric disk. 



3 PHYSICAL REGIME 

3.1 Black Hole Binary Parameters 

The simulations can be scaled, in principle, to any black hole 
mass and orbital separation. In this section, we discuss the 
physical scales for which our simulations could be relevant 
(i.e. physically viable and observationally interesting). 

In general, we are interested in primary MBHs with a 
mass in the range of 10 5 M < M p < 10 9 M . It is not 
clear whether smaller BHs exist in galactic nuclei, and, in 
any case, the radiation from such a low-mass BHs would 
likely be too faint to detect. Likewise, much more massive 
BHs are known to be rare. For the mass ratio q, the interest- 
ing range is 10 -2 < q ^ 1. As we show below, for the set-up 
we study (i.e. with a cavity inserted by hand into the disk), 
the accretion pattern converges as we decrease the mass ra- 
tio to q — 0.01; we expect the results for q < 0.01 would 
remain very similar. In practice, a physical lower limit of 
q ~ 0.01 may also arise from the fact that bound binary 
BHs can be created only in relatively major mergers. In a 
minor merger, the smaller satellite galaxy may not experi- 
ence the torques needed to bring its gaseous nucleus, with 
the low-mass BH, close to the center of the larger galaxy, 
for the BH-BH merger to take place ( |Hopkins et al.| [2006). 
Coupled with the well-established correlations between the 
mass of a MBH and its host galaxy, this suggests that the 
(/-distribution may not extend to values significantly below 
g~0.0lQ 

Following a merger of two MBH-harboring galaxies, the 
MBHs sink to the bottom of the new galactic potential via 
dynamical friction in approximately a galactic dynamical 



timescale (Bcgelma n et al.||1980| . In addition to stellar in- 
teractions (e.g. Preto et al.|201l| , many studies have shown 
that gas in the vicinity of the binary could aid in harden- 



ing the binary down to <C pc separations {e.g. Escala et al 



2005 


|Dotti et al. 


2007| 


Mayer et al.| 


2007 


Lodato et al. 


2009 


|Cuadra et a 


.||2009 


Nixon et al. | 


2011 


Chapon et al. 


2011 


I. We have assumed that such gas in the vicinity of 



the binary cools efficiently and forms a rotationally sup- 
ported thin disk. For a given M v and M„, we are still free 
to choose a physical distance for the orbital radius a, which 
could correspond to a snapshot of the binary anywhere along 
its orbital decay. Accretion disks become self-gravitating, 



5 In principle, less massive BHs may grow from the accretion 
disks around the primary ( jMcKernan et al.|2012^ ; the long-term 
evolution of such systems would be worthy of further study. 
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Table 2. Orbital and residence times for binaries with total mass M = 10 5 M o , M = 1O 7 M , and M = 1O 9 M , at different orbital 
separations a (shown in units of Schwarzs child radius rg ) and mass ratios q = 1 and q = 0.01. The residence time is defined as 



tros = — a(da/dt) _1 , and is adopted from Haiman et al. 2009). The 5 th column approximates the number of orbits until merger 
(N OI b = ires/torbj a * separation a. The last column indicates whether the orbital decay of the binary is driven by gravitational radiation 
or by interacting with the gas disk. Our typical simulation is run for J l 500 orbits, during which the secondary would remain stationary 
to a good approximation. 



M [M ] 


Mass ratio q 


a [rs] 


^orb 


ires [yr] 


N olh 


Migration 


10 5 


1 


10 2 


2.44 hr 


8.0 


2.9 x 10 4 


GW 


10 5 


1 


10 3 


0.46 wk 


2.0 x 10 4 


2.3 x 10 6 


disk 


10 5 


1 


10 4 


14.52 wk 


2.2 x 10 5 


7.9 x 10 5 


disk 


10 5 


0.01 


10 2 


2.44 hr 


2.0 x 10 2 


7.2 x 10 5 


GW 


10 5 


0.01 


10 3 


0.46 wk 


6.2 x 10 3 


7.0 x 10 5 


disk 


10 5 


0.01 


10 4 


14.52 wk 


5.8 X 10 4 


2.1 x 10 5 


disk 


10" 


1 


10 2 


1.45 wk 


8.4 x 10 2 


6.7 x 10 5 


GW 


10 7 


1 


10 3 


0.88 yr 


6.3 x 10 5 


7.1 x 10 5 


disk 


10 7 


1 


10 4 


27.8 yr 


5.7 x 10 6 


2.0 x 10 5 


disk 


10 7 


0.01 


10 2 


1.45 wk 


5.7 x 10 3 


2.0 x 10 5 


disk 


10 7 


0.01 


10 3 


0.88 yr 


1.9 x 10 5 


2.1 x 10 5 


disk 


10 7 


0.01 


10 4 


27.8 yr 


1.5 x 10 6 


5.1 x 10 4 


disk 


10 9 


1 


10 2 


2.78 yr 


8.0 x 10 4 


2.9 x 10 4 


GW 


10 9 


1 


10 3 


88.0 yr 


1.7 x 10 7 


1.9 x 10 5 


disk 


10 9 


1 


10 4 


2.78 X 10 3 yr 


3.0 x 10 8 


1.1 x 10 5 


disk 


10 9 


0.01 


10 2 


2.78 yr 


1.0 x 10 5 


3.6 x 10 4 


disk 


10 9 


0.01 


10 3 


88.0 yr 


9.5 x 10 6 


1.1 x 10 5 


disk 


10 9 


0.01 


10 4 


2.78 X 10 3 yr 


3.0 x 10 s 


1.1 x 10 5 


disk 



and unstable to fragmentation, beyond a radius of order 
> lO 4 (A//lO 7 M )~ 1 r s (where rs is the Schwarzschild ra- 
dius; see, e.g., |Goodman||2003 Haiman et al. 20091. Since 
the binary has to fit inside a gravitationally stable disk, this 
puts an upper limit on the orbital separation. For a fiducial 
binary mass of M = 1O 7 M , we have a max ~ 10 4 rs, or 10 -2 
pc. 

Another practical consideration is the orbital time, 

2tt / M \ ( a \ 3/2 

t orb = ^ = 0.88 ( ) ( ) yr. (8) 



1O 7 M / V10 3 r s 

As we will show, the accretion rate shows periodicity on a 
timescale of ~ t m b- In a realistic survey, it will be feasible to 
look for periodic variations perhaps between 0.1 hr < t OT ^ < 
few yr. Here the lower limit comes from the integration time 
required to measure the flux variations for MBHBs in the 
above mass range (for a survey instrument with a sensitivity 
similar to LSST; Haiman et al.|[2~009 l, and the upper limit 



comes from the duration of proposed time-domain surveys. 
For a M = 1O 7 M binary, this translates to rs < a < 10 3 rs, 
or 10" 6 pc < a < 10" 3 pc. 

For our simulations to be self-consistent, we also require 
a > 100rs, since our Newtonian treatment ignores general 
relativity. Furthermore, at approximately the same binary 
separation, the orbital decay of the binary due to gravita- 
tional wave emission becomes more rapid than the viscous 
time at the edge of the cavity. As a result, the disk de- 
couples from the binary and is 'left behind', rendering our 



jevic & Phinney 



regime (e.g. 


Milosavl- 


|Farris et al.| 


2012| and 



Noble et al.||2012 whose MHD simulations suggest that the 
gas can follow the binary down to smaller separations). 

Finally, throughout our simulations, we fix the bi- 
nary separation; we therefore require that the orbital decay 
should be slow enough for the binary's orbit not to change 
significantly over a few thousand orbits. Assuming that the 



binary is embedded in a thin disk, Hai man et al.| ( [2~009[ ) com- 
pute residence times, t TCS = —a(da/dt)~ , as a function of 
the binary separation a. Among the regimes they consider, 
the one most likely to be relevant for our binaries is the 
so-called 'secondary-dominated' type II migration, i.e. when 
the secondary's mass exceeds the local disk mass. Note that 
there is a trade-off: a longer residence time is desirable since 
it increases the probability of finding such a system; how- 
ever, longer residence times occur at larger separations and 
longer orbital times, which will make it more difficult to 
verify any periodic behavior. 

The residence times and corresponding orbital times, 
for binaries with total mass M = 1O 5 M , M = 10 7 M Q , 
and M = 1O 9 M , and mass ratios q — 1 and q — 0.01, are 
summarized in Table [5] As this table shows, for a 1O 7 M 
binary at separations of a — 10 2 — 10 3 rg, the orbital time 
is between a ~week and a ~year, and the residence time 
is between 10 5-6 years (except for the smallest separation 
of an equal-mass binary, which is already in the GW-driven 
regime, with a residence time of ~ 10 3 years). The binaries 
spend > 2 x 10 5 orbits at these separations, justifying our 
constant a assumption over the < 10 4 orbits we simulate. 
At one extreme of the table, for M = 10 M , q = 1, and 
a = 10 2 rs, the residence time of 8 years is too short for a 
probable detection. At the other extreme, for M = 1O 9 M 
and a = 10 4 rs, the orbital time of 2780 years is much too 
long to pin down over a human lifetime! We also note that 
these results have been derived from simplified migration 
models; recent studies have found significantly longer resi- 
dence times ( Kocsis et al.|2012 b Rafik ov|2012[ ), making the 
above conclusions conservative. 

Finally an important question is whether the binary + 
disk systems simulated here will indeed form a cavity during 
earlier stages of their evolution. For q 0.01, our binaries 
certainly satisfy the standard gap-opening conditions in thin 
disks, requiring that binary torques balance viscous stresses 
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(e.g. Lin & Papaloizou 1986 Artymowicz fc Lubow|[l994| 



On the other hand, the formation of a central cavity, rather 
than an annular gap, also requires that the binary migra- 
tion stalls, allowing the inner disk to drain. Such stalling 
creates the pile-up of gas outside the secondary's orbit, and 
depending on the details of viscosity and binary parameters, 
this pile-up can cause an overflow that fills the cavity. The 



steady-state solutions in Kocsis et al. (2012b see their Fig- 



ures 4 and 5) suggest that for q 0.01, and for a 10 Mq 
binary, a cavity persists down toa~ 300rs , before it is filled 
in; the overflow for smaller [larger] primary BHs occurs ear- 
lier [later]. 

In summary, we conclude that the intuition gained from 
our simulations can be applied to systems with a primary 
BH with M p ~ 10 7 Mq and binary separations lOOrg < 
a < 10 3 rs. For such systems, the binary is in the Newto- 
nian regime, its separation stays constant for > 10 orbits, 
the disk is gravitationally stable, we expect a cavity to be 
present, and the orbital time-scale is conveniently in the pos- 
sibly detectable range (of a week to a year) . Going to more 
massive MBHs, the range of a allowed by requiring the disk 
to be stable is reduced. For M = 10 9 Mq, this leaves a rel- 
atively narrow interesting range around a few hundred rs, 
where the disk is still stable, and the orbital time is a few 
years (see also Tanaka et al.|2012 1 . For less massive MBHs, 
Toomre-stability allows larger separations, out to a ~ 10 6 rs 
for M — 1O 5 M0. The limitations then become long orbital 
times, and too few orbits spent at such large separations. 
Using the models of Haiman et al. (20091, and requiring 



for M p = 10 J M e 



torb < 10 yrs and at least > 10 s orbits 
we find the upper limit a < 10 5 rs (or 10 pc). Addition- 
ally, for such low-mass primary BHs, the cavi ty may close 
at relatively large radii, requiring a > 10 3_4 rs (Kocsis et al. 
2012b I. However as Table [2] shows there remains a sweet- 
spot for 1O 5 M BH's at a ~ 10 4 r s , exhibiting ~ 10 week 
orbital periods and residence times of ~ 10 5 years. 



3.2 Disk Parameters and Accretion Rate 

Our primary goal is to quantify the magnitude and vari- 
ability of accretion across the central cavity. To this end we 
compute the time-dependent accretion rate at the inner edge 
of the simulation domain, r m i n = a 



M(t) = / E(r m i n )v r (r m in)r mill 



(9) 



Since we are unable to track the fate of the gas at smaller 
radii, nor do we allow the masses of the BHs to increase, this 
mass is effectively lost from the simulation. In practice, the 
total mass that is lost is a small fraction of the total initial 
disk mass (at most a few % by the end of each run). 

Because we neglect the self-gravity of the disk, equa- 
tions |T]) are independent of Eo, and there is no unique way 
to assign a physical normalization to M s i m without appeal- 
ing to a disk model which includes additional physics. In- 
stead, MM08 compare the average accretion rate at the edge 
of the integration domain (r m i n = a) to that in a disk in a 
Keplerian potential with the same surface density at r ~ 3a, 



Mbin 

M{ [cc 



M Blm (o) 

67TQ 



H\~ 2 

— J (GMry^-E- 1 ^). (10) 



MM08 find M biu /M ilca ~ 0.2. To make this compari- 
son meaningful, one must assume that Mf ree (3a) ~ M{ lee (a), 
i.e. that the reference, fictitious point-mass disk is in steady- 
state. However, for a steady state disk, specifying a and H/r 
(along with the dominant source of opacity) sets the physi- 
cal value Mfree- For our fiducial values of M = 10 7 Mq and 
a = 10 3 r s (and with a = 0.01 and H/r = 0.1), we find the 
unphysically large value of Mf rcc ~ 10 7 A/Edd (where MEdd 
is the accretion rate that would produce the Eddington lu- 
minosity, with a radiative efficiency of 10%). If we instead 
require Aff rcc ~ MEdd, then this would translate to a much 
thinner disk, with H/r — 4 x 10 -3 . In such a thin/cold disk, 
resolving density waves with the same number of cells as 
in MM08 would require a radial grid resolution ~ 17 times 
higher than in our highest-resolution run, and would be im- 
practical. 

Rather than attempting to compare our binary simu- 
lations to a hypothetical steady-state point-mass disk, we 
chose to leave the surface density normalization So essen- 
tially arbitrary, and instead performed explicit reference 
simulations with a single BH (i.e. q = 0), with the same 
initial conditions as the binary runs. This approach has the 
advantage of explicitly isolating the effect of turning on/off 
the binary. Note, however, that our point-mass runs should 
not be expected produce the steady-state accretion rate for 
a corresponding single-BH system. This is because our ini- 
tial conditions are far from this state, and we do not run the 
simulations long enough (i.e. for a few viscous time at the 
outer edge) to allow it to settle to the correct steady-state. 



4 RESULTS AND DISCUSSION 
4.1 Equal-Mass Binary 

We begin by describing the results of our equal-mass binary 
runs in some detail. Although these are very similar to those 
of MM08, this will serve as a useful point of comparison for 
our unequal-mass runs. 



4-1.1 A Toy Model with Massless Particles 

Before showing the results from our simulations, we consider 
a simple toy model, based on the orbits of non-interacting 
massless test particles around a binary. We populate a 2D 
disk with test particles, centered on the binary's center of 
mass, and leaving out a central cavity. We assign initial 
velocities equal to the Keplerian velocities around a single 
point mass M = M v + M s . We then follow the orbit of each 
test particle in the rotating binary potential, by numerically 
solving the restricted three-body problem for each individ- 
ual particle (for 10 5 particles in practice, using equations 
3.16 and 3.17 in |Murray fc Dermott|2000 l. 

The results of this simple exercise are displayed in Fig- 
ure [2] which shows the locus of the test particles initially, 
as well as after 0.25, 0.5, and 0.75 binary orbits (in a frame 
co-rotating with the binary). As this figure demonstrates, 
there is a tendency for the binary to pull streams of parti- 
cles into the cavity. This, of course, is purely a gravitational 
effect. The toy- model can not be pushed much further in 
time, since after ~ 0.75 binary orbits, the trajectories of 
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Figure 2. The distortion of a disk of massless test particles, 
initially in circular orbits around the center-of-mass of an equal- 
mass binary, with a central cavity. The panels show snapshots of 
the locations of the disk particles after 0.25, 0.5, and 0.75 binary 
orbits, as labeled. The orbits of the particles were followed by 
solving the restricted three-body problem, and are shown in a 
frame co-rotating with the binary. The binary point masses are 
marked by the two [red] dots at x = ±0.5a. The figure illustrates 
the tendency of the binary to create streams of particles entering 
the central cavity, due to gravity alone. 



the test-particles cross - this necessitates a hydrodynam- 
ical treatment. Nevertheless, the figure does suggest that 
an empty cavity can not be maintained by an equal-mass 
binary, even in the absence of pressure or viscosity. Further- 
more, as we will see below, the simulations show accretion 
streams with morphologies quite similar to those in the bot- 
tom right panel of Figure [2] 

4-1-2 Hydrodynamical Evolution: Reaching Steady State 

We next present the results from our equal-mass binary sim- 
ulations. The disk evolves through two distinct stages. As 
explained above, the disk is set up to be out-of-equilibrium, 
and we observe an initial transient state, which lasts for 
~ 1000 orbits. We expect the details of this state to depend 
on the initial conditions. The disk then settles to a quasi- 
steady state, which persists for the rest of the simulation. In 
this quasi-steady state, the disk exhibits significant accre- 
tion, which varies periodically on two timescales, (l/2)tbi n 
and 3tbin- We expect these latter features to be robust and 
insensitive to initial conditions. 

Transient State. Initially, pressure forces and viscous 
stresses act to move material inside of the initial density 
peak at r ~ 10a inward while pressure gradients move ma- 
terial outside of the density peak outward. Once the inner 
disk material reaches r ~ 2a, its surface density becomes 
strongly perturbed and highly non-axisymmetric. Reminis- 
cent of the evolution of the toy model in §4.1.1| two narrow 



mirror-symmetric streams develop, in which the gas flow 
becomes nearly radial. About 10% of the material in these 
streams exits the integration domain at r m i n = a. The rest 
gains angular momentum from the faster moving black holes 
and is flung back toward the bulk of the disk material at 
r ~ 2a. This process maintains a central cavity within the 
disk, with gas streams being pulled in and pushed out on a 
period of ~ 1.5 thia.- It is tempting to associate this period 
with the fibin : fi = 3 : 2 outer Lindblad resonance (LR), 
located at = (3/2) 2//3 a which is the strongest resonance 
available to an equal mass circular binary. However, as also 



emphasized in Shi et al. (20111, the usual interpretation of 



a LR requires the particles at resonance to be on nearly Ke- 
plerian orbits, whereas at ri, gas is flowing in nearly radial 
streams towards the binary. We offer a different geometrical 
explanation for the 1.5 thia variability below. 

As more disk matter flows in from the initial density 
peak to the inner r « 2a region, the streams become more 
dense. When these streams are flung back out and hit the 
opposing cavity wall, they generate noticeable overdensities, 
and deform the circular shape of the cavity to become elon- 
gated and lopsided. These overdensities then rotate at the 
disk's orbital velocity. They spread out and propagate into 
the disk as differential rotation causes them to wind up, cre- 
ating a mirror-symmetric (m = 2), rotating spiral pattern. 

The disk symmetry lends an explanation for the ~ 
1.5 thin accretion rate variability noted above. Fluid at the 
cavity wall (r ~ 2a), and thus the inner edge of the rotating 
spiral pattern, complete an orbit once every ~ 3 thia- Be- 
cause of the m — 2 symmetry, the disk looks the same every 
~ 1.5 thin- Also, the binary returns to the same state every 
1/2 ibin, swapping secondary and primary. This means that 
for the q = 1 binary, the binary disk systems finds itself 
in the same state every ~ 1.5 tbin i-e. streams are generated 
once every 1.5 tbin when the holes come nearest the lopsided 
disk cavity wall. 

Quasi-Steady State. After the initial ~ 1000 orbits, the 
central cavity becomes noticeably lopsided, and stream gen- 
eration becomes stronger on the near side of the cavity. This 
causes more mass to be driven into the cavity wall on the 
far side, pushing that side of the wall still farther away from 
the binary. Since the opposite, weaker stream is generated 
on the now-farther side of the cavity, it is further weakened 
as the disk becomes even more lopsided. This "feedback" 
process continues for a period of 200-300 orbits, after which 
the weaker stream, and its effect on the cavity wall structure 
disappears entirely. The central cavity takes on a steady lop- 
sided shape, with a near-side where streams are pulled from 
the cavity wall by each passage of the holes, and a far side 
where non-accreted material from the streams is flung back 
and crashes into the cavity wall. At the azimuthal locations 
of these crashes, the far-side of the wall develops very strong 
shocks, with Mach numbers up to M w 15. (However, since 
our disk is locally isothermal, this is likely an upper limit 
for the shock strength.) 

We show examples of the two-dimensional surface den- 
sity distributions in Figure [3] The top row of this figure 
shows snapshots at ~ 800 binary orbits (left) and at ~4000 
binary orbits (right), of the inner 6% of the simulated disk 
(i.e. ±6r/a are shown in both directions). The solid cir- 
cle marks the inner boundary of the simulation domain at 
r = r m ; n = a. The larger dotted circle at r ~ 2.08a is 
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q=1.0 q=1.0 




Figure 3. Top row: Surface density distributions for the equal-mass ratio (q = 1.0) binary during a transient, mirror-symmetric state 
after ~ 800 binary orbits (left) and during the quasi-steady asymmetric state after ~ 4000 binary orbits (right). The inset in the top left 
panel zooms in to the inner ±2.5r/a of the disk (at a later orbital phase) in order to show the stream morphology. Bottom two rows: 
snapshots at ~ 4000 binary orbits, during the quasi-steady state phase, for mass ratios q = 0.5,0.1,0.05, and 0.01, as labeled. Each 
panel shows the inner ~ 6% of the simulated disk, extending ±6r/a in both directions. The solid circles mark the inner boundary of the 
simulation at r = r m i n = a. The larger dotted circle at r ~ 2.08a is the position of the (m, I) = (2, 1) outer Lindblad resonance (shown 
only for reference). Surface densities are plotted with the same linear grayscale in each panel, with the darkest regions corresponding to 
a maximum density of 0.4Sq (i.e. the peak of the initial surface density profile). Orbital motion is in the clock- wise direction. 
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the position of the (m,l) — (2, 1) outer Lindblad resonance 
(not present, but shown for reference). The left of these 
two panels illustrates the mirror-symmetric, transient state. 
The zoomed-in inset of this panel shows two weak mirror- 
symmetric streams (at a later orbital phase) reminiscent of 
the streams seen in the toy model of Figure [2] The right 
panel shows the disk after it has settled to its quasi-steady, 
lopsided state, with a single stream. 

The top left panel of Figure|4]shows snapshots of the az- 
imuthally averaged surface density profile of the equal-mass 
binary disk at three different times, after 100, 2000, and 4000 
orbits. For comparison, the density profile is also shown for 
the single point-mass (q = 0) case, after 4000 orbits. As the 
figure shows, the inner circumbinary disk spreads inward 
with time. By 4000 orbits, the disk structure at r/a > 5, 
where the effect of the binary is relatively small, closely fol- 
lows the q = case. However, the density profile remains 
sharply truncated inside r < 2a (i.e. below the point-mass 
case, even at 4000 orbits). The peak density decreases as 
matter drains inward towards the holes and also outward 
towards the outer boundary. The figure also shows that the 
position of the cavity wall moves slightly outward between 
the t=2000 and t=4000 snap-shots. This is because as the 
disk settles to its lopsided quasi-steady state, the cavity size 
grows in the azimuthally averaged sense. 



4-1.3 Torque Balance and the Size of the Central Cavity 

As (the top left panel of) Figure [4] shows, the central cav- 
ity around the equal-mass binary extends to r/a « 2. The 
cavity edge is expected to be located at r co , defined as the 
radius where the negative viscous torque matches the binary 
torque, 



+ 



dT\ 
~dr ) 



= 0. 



(11) 



To gain insight into the transport of angular momentum and 
the clearing of the central cavity, we therefore compute the 
time- and azimuthally-averaged torque densities from the 
binary potential, 

,-2ir 



dT 
dr 



E(r, cj>) — (r, <f>)rd<f> 



bin \"" Jo 
and from viscous stresses, 



dT 
dr 



■ ( '! 

dr 



3 

r v 



(12) 



(13) 



The outer derivative in equation ( 13 1 is taken numerically, 



and all of the above values are measured directly from the 
simulation outputs, except for the binary potential deriva- 



tive in equation ( 12 1 which is given analytically. 



The top row of Figure [5] shows the binary and viscous 
torque densities for q = 1 over the inner 5r/a of the disk, 
during both the initial transient state (left panel) and the 
subsequent quasi-steady state (right panel). There is indeed 
a well-defined central region, where the binary torques ex- 
ceed the viscous torques and can be expected to clear a 



cavity. The transition (computed via eq. Ill is located at 
r cc ~ 1.85 and r — 2.05a in the transient and quasi-steady 
state, respectively, and is marked in both panels by a verti- 
cal dotted line. These vertical lines are also shown in Figure 



[4] and indeed lie very close to radii where the disk surface 
densities remain truncated. 

The small outward drift of the average position of the 
cavity wall from the transient to the quasi-steady state is 
also visible in Figure [4] As the disk transitions to the quasi- 
steady state, the binary torque-density wavelength increases 
in r/a, while keeping approximately the same amplitude. 
These effects can be attributed to the increasingly elongated 
and lopsided shape of the inner cavity; when azimuthally 
averaged, this results in a larger cavity size. 



4-1-4 Accretion Rates 

The most interesting consequence of the lopsided cavity 
shape is on the accretion rate. In the top left pair of pan- 
els in Figure[6] we show the accretion rate, measured across 
the inner boundary of the simulation (r m i n = a) during the 
quasi-steady state, after 4000 binary orbits. The upper panel 
shows the accretion rate as a function of time for « 16 bi- 
nary orbits; the solid horizontal [blue] line shows the time- 
averaged accretion rate during this time, and, for reference, 
the horizontal dashed [green] line shows the accretion rate 
over the same orbits in the q = reference simulation. Inter- 
estingly, the average accretion rate onto the binary exceeds 
that in the single-BH case. The lower panel shows the cor- 
responding Lomb-Scargle periodogram, measured over 100 
binary orbits at a sampling rate of once per simulation time 
step (~ 1000 per orbit). 

As mentioned above, once the quasi-steady state is 
reached, the accretion rate increases from that in the tran- 
sient state by an order of magnitude. As Figure [6] shows, it 
also begins to exhibit strong (factor of ~ 2) variability. Al- 
though this figure samples the accretion only between 4000 
and 4016 orbits, the pattern is remarkably steady, and re- 
peats itself until the end of the simulation. The accretion 
is clearly periodic, and displays two prominent timescales, 
(l/2)ibm and ~ 3it>i n - The stronger variability at one half 
the orbital-time is due to the passage of each black hole 
by the near side of the lopsided disk and the correspond- 
ing stripping of gas streams from the cavity wall. These 
streams are then driven into the opposite side of the cavity 
(as seen in the top right panel in Figure [3] approximately 
135 degrees from the generation point). The second, longer 
timescale corresponds to the orbital period at the cavity 
wall. As mentioned above, when the non-accreted material 
from the streams hits the far-side of the cavity, it creates 
an over-density, which then orbits at the disk's orbital fre- 
quency there of 27r(2.05a) 3/2 (GM) _1/2 « 2.95W Similar 
over-dense lumps have also been found and described in the 



3D MHD simulations of Shi et al. (20111; more recently, 



Roedig et al. (20121 have also mentioned the contributions 



of such lumps to the fluctuations in the accretion rate. 



4.2 Unequal-Mass Binaries 

We next turn to the main new results of this paper, and 
examine the disk behaviour as a function of the mass ratio. 
We start with a qualitative description of how the accretion 
pattern changes as we decrease q, and then present a more 
quantitative interpretation of our results in § |4.3| below. 
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Figure 4. Snapshots of the azimuthally averaged disk surface density at different times (shown by different curves in each panel, from 
100 to 4000 orbits, as labeled), and for different mass ratios (shown in different panels, from q = 1.0 to 0.01, as labeled). In each panel, 
the solid [black] curve shows, for reference, the density profile in the point-mass (q = 0) case after 4000 orbits. The vertical dotted lines 
mark the radius where binary and viscous torques balance (from Figure[5]l; these lie close to the observed cavity edges. The vertical solid 
lines mark the inner edge of the integration domain (r = r m j n = a). In each case, the inner circumbinary disk spreads inward with time, 
but the density profile remains sharply truncated, with a low-density central cavity inside r < 2a. 



4-2.1 Three- Times cole Regime: 0.1 < q < 1. 

The behaviour of systems with mass ratio in the range 
0.1 < q < 1 are illustrated in Figures |3j|4j and |6j These show 
snapshots of the 2D surface density in the quasi-steady state, 
the evolution of the azimuthally averaged density profile, 
and the time-dependent accretion rates, respectively. Addi- 
tionally, in Table [3j and in the corresponding Figure [7] we 
show the time-averaged accretion rate as a function of q. In 
each case, the accretion rate is averaged over 16 orbits in 
the quasi-steady stage, and is quoted in units of the corre- 
sponding rate for a single-BH (q = 0) disk. 

These figures and table illustrate several trends as q is 
lowered from q = 1 -> 0.75 -s> 0.5 0.25 -> 0.2 -s> 0.1: 

(i) The cavity becomes more compact, but less lopsided, 
as one naively expects when the binary torques are reduced. 
These effects are clearly visible in the middle row of Figure|3j 



and also in the corresponding azimuthally averaged density 
profiles in Figure [4] the profiles look remarkably similar for 
q — 0.5 and 0.1, except the cavity for q — 0.1 is smaller, 
and its wall is visibly sharper (as a result of azimuthally 
averaging over a less lopsided 2D distribution). 

(ii) The secondary comes within closer proximity of the 
cavity wall as q is reduced. Again, this is seen in the mid- 
dle row of Figure [3j and occurs both because the cavity 
is smaller, and because the secondary is farther from the 
center-of-mass of the binary. 

(iii) The dense lump, created by the shocks due to the 
"regurgitated" stream-material thrown back out by the bi- 
nary, is still present for q = 0.5 (see the corresponding panel 
in Fig.|3j between 12 and 3 o'clock at r ~ 2a), but is barely 
discernible for q = 0.1. Again, this trend is unsurprising - as 
the torques diminish, one expects weaker shocks and smaller 
overdensities in any resulting lump. 
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Figure 5. Azimuthally- and time-averaged torque density profiles in the inner disk for the equal— mass binary (top two panels) and for 
unequal-mass binaries (other panels, with different mass ratios q as labeled). The top left panel corresponds to the mirror-symmetric 
transient stage (after ~ 800 orbits) and the top right panel to the asymmetric quasi-steady stage (after ~ 4000 orbits). Only the quasi- 
steady stage is shown for the q < 1 cases. In each panel, the dashed [black] curves show the gravitational torques from the binary, and the 
red [dot-dashed] curves show the negative viscous torques. The vertical dotted line marks the radius where the viscous and gravitational 
torques balance (equation |11| ; these are close to where the azimuthally-averaged surface density profiles are found to be truncated (see 
Fig.Ek. Time averages are taken over 25 orbits at a sample rate of 20 per orbit. 
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Figure 6. The time variable accretion rate across the inner boundary of the simulation measured at r = r m i n = a (top of each pair of 
panels) and the corresponding Lomb-Scargle periodogram (bottom of each pair). Panels are displayed in order of decreasing binary mass 
ratio, starting from q = 1.0 at top left (on the previous page) to q = 0.01 at the bottom right (on this page). The average accretion rate in 
each panel is denoted by the solid [blue] horizontal line. The dashed [green] horizontal line in each plot shows the average accretion rate 
for the point mass (q = 0) case for reference. For q > 0.05, the accretion rate is strongly modulated by the binary, with either one, two, 
or three distinct periods present simultaneously, depending on the value of q (see text for detailed explanations). For q < 0.05, the binary 
still reduces the mean accretion rate noticeably, but does not imprint strong periodicity; the q = 0.01 binary is nearly indistinguishable 
from a single BH. 



Table 3. The mean accretion rate Mbim averaged over 16 orbits in the quasi-steady state, for binaries with different mass ratios. The 
rates are shown in units of the corresponding rate M„q found in a single-BH (q = 0) simulation. The four rows show results for different 
combinations of radial and azimuthal resolutions. The second row is our fiducial resolution. 
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(iv) The accretion streams in the q = 0.1 panel in Figurej3] 
become noticeably weaker. However, the average accretion 
rate, shown in Figures [6] and [7] and in Table [3] stays at 
approximately ~ lAM q o (i.e. exceeding the single-BH ref- 
erence value) over the range 0.5 < q ^ 1. For q < 0.5, the 
average accretion rate drops more rapidly, falling by nearly 
a factor of two to ~ 0.76sM 9 o by q = 0.1. 

(v) Despite a drop in the average accretion rate, the max- 
imum accretion rate (i.e., the amplitude of the spikes in 
Fig. [6| increases as the mass ratio is decreased to q ~ 0.1. 
For example, at q = 0.1, the fluctuations in the accretion 
rate reach a factor of ~ 4, compared to a factor of ~ 2 in 
the q = 1 case. As a result, even though the average accre- 
tion rate is ~ 50% below the q = 1 case, the largest spikes 



reach ~ 25% above the corresponding maxima for q = 1. 
This could make it easier to detect and identify such bina- 
ries, despite their weaker average accretion. 

(vi) Perhaps the most interesting result is shown by the 
Lomb-Scargle periodograms in Figure [6] As q decreases, 
power is clearly traded from both the (l/2)tbi n and the 3tbin 
variability timescales into the tbin timescale. This is because 
of the increased proximity between the secondary and the 
cavity wall, and a corresponding larger distance between the 
primary and the cavity wall noted above. As a result, as q is 
decreased, the secondary begins to dominate the variability, 
pulling stronger accretion streams off the cavity wall once 
per binary orbit. 
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Figure 7. The average mass accretion rate as a function of q, 
and for simulations with different spatial resolution. This is the 
same information as in Table [3] except here shown graphically. 

Focusing on the last finding: in the 0.1 < q < 1 case, we 
find that the time-dependent accretion rate displays three 
distinct and sharply defined periods, with well-defined ratios 
at 0.5, 1, and 3£bin- While the last of these reflects the orbit of 
the dense lump at r ~ 2a and could depend on details of the 
disk properties, the first two periods are fixed by the binary 
alone and are independent of the disk. The 1:2 period ratio 
is therefore a robust prediction; if observed, it could serve as 
a smoking gun signature of a binary. 

4-2.2 Single-Orbital-Timescale Regime: 0.05 < q < 0.1. 

As q is further decreased, the overall distortions to the disk 
become less pronounced, and approach a nearly axisymmet- 
ric, tightly wound spiral pattern (see the q — 0.05 panel 
in Figure [3]| . The distance between the cavity wall and the 
secondary further shortens, and the accretion variability be- 
comes dominated entirely by the streams created by the 
secondary's passage once per orbit. As Figure [5] shows, for 
q — 0.075 and 0.05, the accretion rate displays a nearly 
perfect sinusoidal variation, with the corresponding Lomb- 
Scargle periodograms showing a single spike at the orbital 
timescale of the binary. For q = 0.075, the fluctuations are 
still large (factor of ~ 3), but by q = 0.05, the fluctuations 
strongly diminish. Interestingly, however, we find that the 
average accretion rate is most strongly suppressed among 
all of our runs for q = 0.05 (by a factor of ~ 5 compared to 
the single-BH case). 

As mentioned above, the disappearance of the (l/2)tbm 
variability timescale is easy to understand qualitatively: 
once the primary BH remains very close to the center-of- 
mass, its compact orbital motion no longer impacts the disk 
far away. The disappearance of the 3ibi n variability timescale 
is also clearly attributable to the lack of any dense lump near 
the cavity wall for q < 0.1. However, the reason this lump 
disappears is less obvious, and warrants some discussion. 

Figure [6] shows that as q decreases, the power in the 
(l/2)tbin and 3tbin timescales decreases in tandem, main- 



taining a roughly constant ratio. This suggests that the vari- 
ability on these two timescales are causally linked. Such a 
causal link could arise either through the strength or the 
frequency of the accretion streams: 

(i) As q decreases, the cavity becomes less lopsided, and 
the accretion rate spikes become weaker. This suggests that 
when these weaker accretion streams are flung back to the 
cavity wall, they create less overdense lumps. For q ^ 0.1, 
the lump may not survive shear stresses and pressure forces, 
and may dissolve in less than an orbital time. 

(ii) As can be eyeballed from Figure [3] the stream impact 
zone (dense region outside the cavity) extends over an az- 
imuth of A<f) ~ 100 — 120° . The orbital period at the cavity 
edge is ~ 3t^ n , implying that multiple streams can hit parts 
of the same lump if and only if streams are generated more 
than once per binary orbit. 

To test whether the strength or the frequency of the 
streams is more important for lump generation, we repeated 
our simulation for an equal-mass binary, but we placed one of 
the BHs artificially at what would be the real binary's center 
of mass. The second hole still orbits at r = a/2 as usual. In 
this setup, the cavity wall is perturbed by streams with a 
similar strength as in the real q — 1 simulation, but now only 
once, rather than twice per orbit. We found that this "one- 
armed" binary did not create a significant over-density at 
the cavity edge, nor did it excite a significant elongation of 
the cavity. We therefore conclude that multiple, overlapping 
streams are required to generate a strong lump that survives 
for an orbital time. This explains the roughly constant ratio 
of power between the (l/2)tbin and 3tbin timescales: as the 
power at (l/2)tbm is reduced, the lumps become weaker, 
decreasing the power at 3£bin- This result also explains why 
the cavity becomes less lopsided as q decreases. 

4-2.3 Steady-Accretion Regime: q < 0.05. 

As we continue to decrease q from 0.05 to 0.01, we find 
yet another distinct regime. The overall morphology of the 
snapshot of a q — 0.01 disk in Figure [3] looks similar to the 
q — 0.05 case, except the nearly-concentric perturbations 
are even weaker, and the cavity still smaller. However, the 
similarity is quite deceptive, with the movie versions of these 
figures showing a striking difference]^] In the q = 0.05 case, 
accretion streams form and disappear periodically, but in the 
q = 0.01 case, the disk pattern becomes constant and un- 
changing (in the frame co-rotating with the binary). There 
is still a visible accretion stream, hitting the inner boundary 
of the simulation just ahead of the secondary's orbit, but the 
stream steadily co-rotates with the binary. 

As Figure [6] shows (see also Fig. [7] and Table [3J, the 
average accretion rate has reached its minimum at q = 0.05. 
As q is further decreased, the accretion rate becomes steady, 
with no fluctuations, and its value increases back towards 
the q = rate (dashed horizontal [green] line in Figure [6|. 
Essentially, for such a light secondary, the system begins 
to resemble a disk with a single BH. Although a cavity is 
still clearly present (moving from r ~ 1.6a at q = 0.05 to 

6 Movie versions of the snapshots in Figure [3] are available at 
http://www.astro.columbia.edu/~dorazio/moviespage 
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r ~ 1.4a at q = 0.01), it is being refilled, as for a single BH. 
In fact, after 4000 orbits, the azimuthally averaged density 
profiles are nearly identical for g = and q = 0.01 (see 
Fig. [4}. 

Although the secondary still excites small but visible 
ripples in the disk (Fig. [jjl, by q = 0.01 it can no longer 
exert a large enough torque to pull in large streams and drive 
them back out to produce a lopsidedness in the circumbinary 
disk. The ripples are deeply in the linear regime, and they 
resemble the tightly-wound spiral density waves launched 



in protoplanetary disks (e.g. Goldreich & Tremaine 1980 
Dong et al.||2011| |Duffell fc MacFadyen||2012[ ), except here 



our background disks have a pre-imposed central cavity. In 
the deeply linear regime, the waves are excited by resonant 
interactions with the disk, and non-linear coupling to an 
m — 1 mode is no longer possible. Note the disappearance 
of any ripples in the azimuthally averaged surface densities 
for q < 0.05 (Fig. [If. 

4.3 Summary and Discussion 

In summary, we find that the behavior of the accretion rate 
across the circum-binary cavity as a function of q can be 
categorized into four distinct regimes: 

(i) Two-timescale regime; q = 1. Confirming previous re- 
sults, an equal-mass binary maintains a central low-density 
cavity of size r ~ 2a. However, the time-averaged accretion 
rate exceeds that for a point-mass case; there are factor of 
~ 2 fluctuations around the average on two prominent time- 
scales, (l/2)t b in and « 3tbin- 

(ii) Three-times cale regime; 0.1 < q < 1. The time- 
averaged accretion rate drops by a factor of ~two by q = 0.1; 
however, the fluctuations continue to grow and the peak 
accretion rates remain high. There are three time-scales 
present, (l/2)t b in, t hin , and « 3t b i n . 

(iii) Single-orbital-timescale regime; 0.05 < q < 0.1. In 
this regime, the average accretion rate and its fluctuations 
both drop with decreasing q. Variability is dominated by 
the secondary alone, and is nearly sinusoidal on the binary 
period ibm- 

(iv) Steady-Accretion regime; q < 0.05. The accretion be- 
comes steady, and the system overall resembles a cavity- 
filling single-BH disk, with small perturbations due to the 
secondary in the linear regime. 

4-3.1 Understanding the q-Dependence of the Disk 

To better understand the reason for the transitions through 
each regime, we note that lowering the binary mass ratio has 
two main consequences. First, the position of the secondary 
(primary) moves away from (towards) the binary's center of 
mass, 



r s (q) = o(l + q)~ 



r p (q) = a(l + l/q)~ 



(14) 



Second, the size of the central cavity decreases. This 
is apparent visually from both the 2D surface density snap- 
shots (Fig.[3| and the azimuthally averaged profiles (Fig.[4|. 
The latter figure also shows the locations of the cavity edge 
r co expected from balancing the azimuthally averaged grav- 
itational and viscous torques (Fig. [5J , which agree well with 
the observed cavity sizes. In the top panel of Figure [§] we 



explicitly show r ce as a function of q. The points in the figure 
have been obtained by balancing the azimuthally averaged 
viscous and gravitational torques measured in the simula- 



tion (equation 111. The black line is an empirical fit to the 



(15) 



data points at fiducial resolution, given by 

r cc (q)~ 3.0 In (q 01 + l) a. 

Figure |8j shows, in particular, that for mass ratios 0.5 < 
q ^ 1.0, the average cavity size stays approximately con- 
stant, while for mass ratios q < 0.5, the cavity shrinks more 
rapidly with q (cf. |Artymowicz fe Lu bow 1994; Milos avTjevic| 
fc Phinney|2005 I. We attribute this behavior to weaker out- 
ward binary torques for smaller q. The binary torque, as a 
function of q and disk coordinates, can be computed from 
the analytic binary potential (eq. [2]| and is given by 



Tbin (r, (p, q) 



-GM 



+ 



a 

1 + q- 1 



ar sin (cf> — fibin^) 
(1 + q) (1 + q- 1 ) 5 

3/2 



2ar cos (<f> — ilbini) 
1 + q- 1 



2 , 

r + 



l + < 



2ar cos (0 — fibini) 
1 + ? 



-3/2' 



(16) 



We plot Tbin(r, <j>, q) in the middle panel of Figure [8] at 
three constant radii in the disk, at an azimuthal coordinate 
~ 26° behind the position of the secondary, where the out- 
ward torque is maximal. The figure shows that the outward 
torque, for example at r = 2a, increases weakly from q = 1.0 
to 0.5, but then decreases more rapidly, below a l/3rd of its 
peak value, as q decreases to q = 0.05. Thus, the cavity size 
as a function of q is well correlated with the strength of the 
outward binary torque which clears the cavity: both are ap- 
proximately constant or weakly increasing as q drops from 
1.0 to 0.5, and then both decrease rapidly as q falls towards 
0. 

We can next use the apparent cavity size r ce (q) and the 
secondary's location r s (q) to estimate the impact of each 
hole as it passes by a fluid element in the neighboring cavity 
wall. We first evaluate the radial velocity imparted to the gas 
by the passage of the secondary, invoking the impulse ap- 
proximation. This assumes that the azimuthal deceleration 
just before and acceleration just after the passage cancel, 
and only a net radial velocity component is imparted onto 
the fluid element. This radial velocity is approximated as the 
product of the radial force from the BH at closest approach, 
and the duration of the encounter. Here we ignore the non- 
axisymmetric elongation of the cavity, and take the impact 
parameter at the closest approach to be the distance between 
the secondary and the average cavity wall edge r cc — r a . We 
note that r ce is a good approximation for the distance from 
binary center of mass to the near side of the cavity in each 
of our runs. This can be eyeballed by comparing the values 
of r ce in Figure [8] to the 2D density profiles in Figure [3] We 
assume that the corresponding encounter time is given by 



^bin^s f^ce^c 



(17) 



where the denominator is the relative velocity of the disk 
material and the passing secondary. During this encounter, 
the primary is passing by in the opposite direction, at a 
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Figure 8. Cavity size and binary torques as a function of mass 
ratio q. The top panel shows the position of the cavity wall as 
a function of q, in runs with different radial and azimuthal res- 
olutions, as labeled. The points mark the radii r cc at which the 
azimuthally averaged viscous and gravitational torques balance 
(equation The black line is an empirical fit to these data 

points. The middle panel shows the binary torque (eq. |16| l eval- 
uated at three constant radii and at an azimuth ~ 26° behind 
the secondary, where outward torque is maximal. For q < 0.4, the 
binary torques decrease steeply with lower q. The bottom panel, 
shows an estimate of the radial velocity imparted to a disk ele- 
ment at the cavity wall when the secondary (blue line) and the 
primary (black line) pass by. The estimate is based on the empir- 
ical fit in the top panel, and on the impulse approximation. 



distance r cc + r p , with a corresponding encounter time of 

At ^ n h r rtn r ■ (18) 

" ^ bin ' p ~ Ji cc' cc 

The net increase in radial velocity due to the secondary pas- 
sage on the near side is then 



5v s = 



GM(l + g- 1 )- 1 At s 



+ 



GM{l + q y l l\t p 



[r ce (q) - a(l + q)- 1 ]^ [r cc (q) + a(l + q- 1 



(19) 



while the net increase in radial velocity due to the pri- 
mary passage Sv p is found by swapping p's and s's in equa- 
tions JT7l) and fl8l). 



In the bottom panel of Figure |8j we show 6v B and Sv p 
as a function of q. First we see that as q is decreased from 
1, the impact of the primary passage by the near cavity wall 
decreases monotonically. This is simply due to the increase 
in r cc — r p as q decreases, which lessens not only the radial 
force felt at the cavity wall but also the encounter time. 

Next, notice that as q is decreased from 1 to ~ 0.4, the 
impact of the secondary passage increases considerably. This 
is because, opposite to the primary case, r ce — r s decreases 
with q. The most important effect of this is to increase the 
encounter time. Even though r cc — r a decreases, the sec- 
ondary's mass decreases more quickly and the radial force 
exerted by the secondary lessens with decreasing q. However, 
as the secondary moves closer to the cavity edge with de- 
creasing q, its orbital velocity increases as r 2 oc l/(l+q) 2 . At 
the same time the cavity edge moves inwards and the orbital 
velocity of fluid elements there also increases but at a slower 



rate proportional to r co oc [31n (q° 1 + l)] H *. The re- 
sult is a matching of the two orbital velocities at q ~ 0.4 
and a formally infinite encounter time at that mass ratio. 
In reality, however 5v s will be cut off as a stream is gen- 
erated or as the fluid element moves to larger r along the 
non-circular cavity edge. If we cap the encounter time at 
~ l/2tj 3m (roughly a stream generation time), then the very 
large Sv s near q ~ 0.4 is truncated. With this restriction, 8v s 
decreases for all q > 0.1 and the secondary reaches its largest 
perturbing power at q ~ 0.1. As q decreases below this turn 
around mass ratio, the impact of the secondary passage de- 
creases as the encounter time again becomes shorter and the 
radial force due to the secondary continues to decrease. 

The above analysis elucidates our results in the three- 
timescale regime and the transition to the single-timescale 
regime. As q decreases from 1 to 0.1, less disk fluid is pulled 
into the cavity over the course of a whole binary orbit, caus- 
ing a drop in the time-averaged accretion rate. However, due 
to its closer proximity to the cavity wall and a longer en- 
counter time with fluid elements at the cavity wall, the sec- 
ondary becomes more effective at generating brief spikes in 
accretion (while at the same time, the effect of the primary 
diminishes). As a result, the maximum accretion rate, during 
the narrow spikes that occur once per binary orbit, increases. 
This latter trend changes for q < 0.2, where the encounter 
time between a disk fluid element and the secondary quickly 
decreases and the secondary becomes too light to generate 
strong accretion spikes. At around the same q ~ 0.2, the 
primary's orbit becomes too compact to generate noticeable 
time-variable perturbations at the cavity wall. 

The impulse analysis illuminates the effect that q has 
on the magnitude of the induced accretion streams as well 
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Figure 9. The binary torques (from eq. |16| , evaluated at the 
cavity wall edge (r co (q); from eq. |15| l. The top panel shows con- 
tours of the total (primary+secondary) torque as a function of q 
and azimuthal angle <j>. The secondary is located at <j> = 7r. Out- 
ward [inward] torques are shown in red [blue] shades. The bottom 
panel shows a one-dimensional slice of the torque per unit disk 
angular momentum at an azimuth indicated by the horizontal 
line in the top panel, just in front of the secondary (A<f> ~ 26°). 
These inward torques are largest at a mass ratio of q = 0.2. This 
suggests that the strongest streams can be pulled into the cavity 
(and flung back out) near this binary mass ratio, as seen in our 
simulations. 



as the relative importance of the secondary and primary for 
generating streams. However, the radial velocities are on the 
order of the relative velocities between the perturbing hole 
and the disk fluid element (hence stream generation), mak- 
ing the impulse approximation suspect. For this reason, we 
next investigate the azimuthal component of the accelera- 
tion imparted by the binary, by computing the torques from 



gle (p. Red regions correspond to outward torques and blue 
regions corresponds to inward torques. The secondary is po- 
sitioned at <f> — 7t and moving in the direction of increasing 
<f>. In the bottom panel of Figure [5] we show the torque per 
unit disk angular momentum at an azimuth just in front of 
the secondary (A<f) ~ 26°, where inward torques are largest). 
These inward torques are largest at a mass ratio of q ~ 0.2. 
This suggests that the strongest streams can be pulled into 
the cavity (and flung back out) near this binary mass ratio. 
This is consistent with our findings that there is a dimin- 
ished (relative to either q = 1 or q = 0) time-averaged accre- 
tion rate, but a maximum in the peak accretion rate, near 
q ~ 0.2Q More generally, this is further evidence that the 
general trend in accretion magnitudes can be understood 
simply via the kinematic interaction of the binary potential 
with the disk, given the geometry and size of the cavity. 

Figure [3] shows that the degree of elongation of the cav- 
ity decreases with mass ratio for q < 0.5. As we showed in 
§4.2.2| this is at least partly due to the decrease in stream 
overlap at impact, because of the diminishing strength of 
streams generated by the primary. However, the middle 
panel of Figure [8] shows that the outward torques, at po- 
sitions just behind the secondary have a maximum near 
q ~ 0.4. Such outward torques are responsible for direct- 
ing the streams outward, and the strength of their impact 
into the cavity wall. We therefore expect that the decrease in 
these outward torques below q < 0.4 also contributes to the 
diminishing lopsidedness of the cavity and the overdensity 
in the lumps that appear at the impact sites. 



4-3.2 Viscosity Dependence and Comparison with MM08 

Our main qualitative conclusions for the equal-mass binary 
case, including the morphology of the disk, and the accretion 
rate, are in good agreement with MM08. Nevertheless, we 
find two small discrepancies with MM08, which are related 
to different treatments of viscosity. 

First, as can be seen by comparing MM08's Figure 1 
with the top left panel of our Figure [4j the azimuthally av- 
eraged density of the MM08 disk near r ~ 2a increases with 
time as the fluid builds up against the cavity wall. Roughly 
speaking, MM08 find that the peak surface density barely 
changes from its value imposed in the initial condition; only 
the location of the peak drifts inward, as the inner regions 
of the disk are accreted toward the binary, but then held 
up near r ~ 2a. In contrast, we find that the peak surface 
density diminishes over time; i.e. our disk drains too rapidly 
too allow a significant build-up in the surface density near 
r ~ 2a. 

Second, comparing MM08's Figures 7 and 8 with the 
top left panel of our Figure [6] we see that although the 



equations (15 1 and (16 1 



The top panel of Figure [9] shows contours of the total 
torque, summed over both binary components, evaluated at 
the cavity wall edge, as a function of q and azimuthal an- 



7 |Hayasaki et al.||2007[ l find that the accretion rate onto the holes 
increases by a factor of ~ 2 (a factor of 3.25 for the secondary and 
no change for the primary) for a circular binary with mass ratio 
q = 0.5 compared to an equal mass binary. However, their simu- 
lations were run for only for 40 orbital times. Since the behaviour 
that leads to the observed accretion rates here does not arise un- 
til after ~ 1000 orbital times (and also since the accretion rate 
depends on gas dynamics within the region near r ~ 1.8a, where 
Hayasaki ct al. (2007) inject mass as a boundary condition), it is 
difficult to make a comparison between the two studies. 
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Figure 10. Accretion rates at the inner boundary r m i n = a 
for the equal-mass binary as in the top left panels of Figure [6] 
except for a viscosity coefficient a reduced by a factor of two, from 
0.01 to 0.005. The solid horizontal [blue] line in the top panel 
is the average accretion rate, and the bottom panel shows the 
Lomb-Scargle periodogram over a duration 3 times longer than 
displayed in the top panel (as in MM08). Comparing to Figure[6] 
we see that a stronger viscosity leads to a more regular accretion 
pattern. 



magnitudes of the time-averaged accretion rates are simi- 
lar, the detailed variability and the resulting Lomb-Scargle 
periodograms are not identical. 

To investigate the cause of these discrepancies, we re- 
ran our q — 1 simulation with the value of the viscosity 
parameter a reduced by a factor of two, from 0.01 to 0.005. 
The accretion rate found in this case is shown in Figure fl0| 
Comparing with the top left panels of Figure [6] we see that 
a weaker viscosity leads to a less regular accretion pattern, 
with the period at (l/2)£bin less dominant, and a handful of 
new frequencies present (with an especially strong new peak 
appearing at ~ 5tbin). Comparing our lower-a results with 
MM08's Figures 7 and 8, we find that the fluctuations in 
the accretion rate and the Lomb-Scargle periodogram match 
very closely those of MM08. However, in this case, the overall 
magnitude of the accretion is ~ 2 — 3 times lower than in 
MM08. This can be attributed to our factor of ~ 2 — 3 lower 
density near the cavity edge (which remains for a = 0.005). 

The appearance of new frequencies in the periodogram 
of Figure |10| compared to the top left panel of Figure |6j 



can be explained by considering the Fourier-transform (FT) 
of each accretion time series. Imagine a train of accretion 
spikes similar to those in the bottom left of Figure [6] Such 
a time series can be constructed by convolving the shape 
of the individual accretion spikes with an impulse-train of 
Dirac-delta-functions. By the convolution theorem, the FT 
of this convolution is the product of the FT of the impulse 
train and the FT of the spike shape. The FT of an impulse 
train with period T is another impulse train in the frequency 
domain spaced by T~ l . The FT of the spike shape will be 
qualitatively similar to a sine function which tapers off and 
rebounds as it moves away from the fundamental frequency. 
The width of this sine function will vary inversely with the 
temporal width of the spike. Thus our periodogram should 
resemble a train of Dirac-delta-functions spaced by T , but 
with amplitudes modulated by a sine function. The q = 1, 
a — 0.01 accretion time series can be considered as a product 
of two such impulse trains, one at Ti/ 2 = (l/2)tbin, and one 
at T3 =~ 3tbin- Likewise, the q = 1, a = 0.005 case is a 
product of a = (l/2)ibin, and a T5 =~ 5£bin train. The 
differences in the two periodograms can then be attributed 
to the different shapes of the time domain envelopes of the T3 
vs T5 periods. The shape of the T5 period spike (in the a = 
0.005 case) is narrower and thus the sine function envelope 
in the frequency domain is wider, allowing us to see the 
'extra' frequency spikes, spaced by T 5 _1 . 

In summary, we have found good overall agreement with 
MM08, but were unable to reproduce a perfect match to 
their time-dependent accretion rate for a single choice of a: 
we can either reproduce their periodogram (for a = 0.005), 
or their time-averaged accretion rate (for a = 0.01). We 
suspect that these discrepancies are caused partly by dif- 
ferences between FLASH2 (used in MM08) and FLASH3 
(used here). In particular, the two versions implement two 
different routines to compute viscous fluxes. The FLASH3 
viscosity routine lias been updated by the method detailed 
in (Edgar 2006) to conserve angular momentum. More im- 
portantly, we also conclude that the details of the accretion 
rate variability depend on the magnitude of disk viscosity 
- at least for a < 0.01. This suggests that once the two 
periods at t OI b and (l/2)i or b are robustly identified in an 
observed system, the rest of the features in the periodogram 
can probe disk properties, such as the viscosity. We save the 
full exploration of viscous effects for a future study and here 
continue to run all other simulations at the fiducial value of 
a = 0.01. 



4-3.3 Resolution Study 

Up to now we have discussed the results from our fidu- 
cial set of medium-radial, low-azimuthal resolution runs 
( [Mid Ar, Lo A(j)] in Table [y. Ideally, we would repeat these 
runs at increasingly high radial and azimuthal resolutions, 
until the results converge. Unfortunately, this is computa- 
tionally prohibitive, and we instead chose the following ap- 
proach. 

(i) For the q = 0.1 case only, we performed two higher- 
resolution runs ([MidAr, MidA0] and [HiAr, HiA^] in Table 
[T| and one lower resolution run ([LoAr, LoA0] in TableJT]) to 
look for signs of convergence. The q = 0.1 case is at the cusp 
of the three-timescale and single-timescale regimes, where 
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the accretion behaviour is particularly sensitive to q. We 
expect it to be the most sensitive to resolution effects, as 
well, and therefore chose this case for our resolution study. 

(ii) We repeated each of our runs at the lowest resolution. 
Redoing the entire set of runs allows us to assess whether dif- 
ferent g's are affected by resolution differently. Furthermore, 
at the faster speeds at lower resolution, we could afford mul- 
tiple and longer runs at the cusp between the variable and 
non- variable regimes in order to determine the nature of this 
transition. 

Figure [Tl] gives a visual impression of the surface den- 
sity distribution at the four different combinations of reso- 
lutions for q — 0.1. While the images look similar overall, 
they show a clear trend: as the resolution is increased, the 
lumps near the cavity wall become sharper and more over- 
dense, and the cavity becomes larger and more lopsided^ 
This trend can be attributed to numerical dissipation when 
the regurgitated streams impact the cavity wall. Increasing 
the resolution implies weaker numerical diffusion, stronger 
accretion streams, more momentum carried into the disk, 
and an overall more efficient driving of the m = 1 mode. 

Figure [12] shows the corresponding time-dependent ac- 
cretion rates at the inner boundary of the simulation. 
The accretion pattern look visibly different for the lowest 
resolution run which exhibits a maximum accretion rate 
~ 3 times lower than its higher resolution counterparts. 
However, encouragingly, the difference between our fidu- 
cial [MidAr,LoA<^] and the highest-resolution [HiAr, HiA</>] 
cases is modest. The differences in the Lomb-Scargle peri- 
odograms can again be explained by considering the Fourier- 
transform of each accretion time series as discussed above 
in j pT3~2l 

Tableland the corresponding Figure [7] shows the time- 
averaged accretion rates as a function of q at different res- 
olutions. In all cases, at the same fixed q, we find that in- 
creasing resolution produces a higher accretion rate. This is 
consistent with the interpretation above that better resolu- 
tion allows stronger accretion streams. Interestingly, we find 
a strong correlation between the values of Mhin/M q o listed 
in Tableland the accretion patterns seen in Figure [12] runs 
at different resolutions but with similar values of Mbm/Af 9 o 
have very similar accretion patterns (including the variabil- 
ity and the values of the maxima). The result of increasing 
[decreasing] the resolution can therefore be interpreted as a 
shift of the accretion behavior to lower [higher] mass ratios. 

Encouragingly, the two higher azimuthal resolution runs 
at q = 0.1, also plotted in Figure [7] lie closer to each other 
than the two lower resolution runs. Since the resolution steps 
are evenly spaced, we consider this evidence that the simu- 
lations are converging monotonically with resolution. 

To find whether there is a cut-off mass ratio between 
variable and non-variable accretion (or whether our simula- 
tions do not run long enough to capture variable behavior 
for small q) we have run a number of low-resolution simu- 
lations spaced by Aq — 0.005 (see Table At our fiducial 
value of a we find the transition from periodically variable 
to steady accretion occurs at q ~ 0.06. At the mass ratio of 



8 Figure [S] shows that the average position of the cavity wall, 
found from the azimuthally averaged torques in equation < |11[ |, is 
largely unaffected over the range of resolutions studied here. 



q = 0.065, the disk becomes eccentric, and causes variable 
accretion after only 3000 binary orbits. The q = 0.06 case is 
run for 8000 orbits and never becomes eccentric. 



4-3.4 Caveats 

However instructive, the simulations presented here are still 
of course simplified models of a real binary disk system, and 
it is worth listing some major caveats. 

(i) Our simulations are two dimensional - we expect that 
the 3D vertical structure could modify the structure of the 
accretion streams, including their g-dependence. 

(ii) Our simple a-viscosity prescription may be inaccu- 
rate, especially in the nearly radial accretion streams. Re- 
cent MHD simulations in the q = 1 case find a much larger 



effective a than we adopted here ( Shi et al. 2011 Noble et al 



|2012p . We have argued that much of the accretion behavior 
can be understood by gravitational kinematics alone, but 
this needs to be verified explicitly, especially in light of the 
a-dependence of our results. 

(iii) We assumed a locally isothermal equation of state; 
more realistic equations of state could have an especially 
large impact on the strength and dissipation of shocks at 
the cavity wall. 

(iv) Our initial conditions correspond to an unperturbed, 
near-Keplerian, circular disk, and circular binary orbit, with 
a significant pile-up of gas. In reality, accretion onto the bi- 
nary could produce significant binary (as well as disk) eccen- 
tricity (e.g. jCuadra et al.||2009| |Lodato et aT1|2009} |Roedig| 



et al. 20111. In this case, the variability we found would 



most likely have a more complex structure (e.g. Hayasaki 
|et al.[2 007) due to the plethora of resonances available in an 



eccentric binary potential (e.g. Artymowicz &: Lubow|1994 | 

(v) Although we begin with an "empty" cavity, this cavity 
may overflow already at a large radius ( Kocsis et al.|2012b I . 
Future studies should construct a self-consistent initial den- 
sity profile, by evolving the binary's orbit from large radius 
through gap clearing, and thus determining whether a true 
pile-up occurs. 

(vi) We have ignored the radiation from the gas accreted 
onto the BHs. Given that we find high accretion rates - 
comparable to those for a single BH - the secondary BH 
can be fed at super-Eddington rates, and the flow dynamics 
can then be strongly affected by the radiation. 

(vii) We have not allowed accretion onto the BHs, and 
have excised the inner region r < a from the simulation do- 
main. This could have an impact on the dynamics of the 
streams that are flung back towards the cavity wall, and 
therefore on the formation of the dense lump, the lopsided- 
ness of the cavity, and the variable accretion patterns. 

These caveats should all be pursued in future work, to 
assess the robustness of our results. We expect that our main 
conclusions, namely that the accretion rate is strongly mod- 
ulated by the binary, and that the power-spectrum of the 
accretion shows distinct periods, corresponding to the or- 
bital periods of the binary and the gas near the location of 
the cavity wall, will be robust to all of the above caveats. 
However, the numerical values, such as the mean accretion 
rate, and the critical value of q for the transition between 
variable and steady accretion, will likely be affected. 
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q = 0.1 [MidAr, LoA<f>] 




q = 0.1 [MidAr, MidA^] 



q=0.1 [HiAr, HiA0] 




Figure 11. Two-dimensional surface density distributions at 4000 orbits, as in Figure[3] except for the single mass ratio q = 0.1, and for 
four different combinations of high/low radial and azimuthal resolutions, as labeled. Increasing the spatial resolution decreases numerical 
diffusion, leads to sharper features, and allows stronger accretion streams. The stronger streams lead to more overdense lumps where the 
regurgitated streams hit the cavity wall. As a result, the cavity becomes larger and more lopsided as the resolution is increased. 



5 CONCLUSIONS 

We have investigated the response of an accretion disk to 
an enclosed binary via two-dimensional, Newtonian, hydro- 



dynamical simulations. As previous work has shown (Arty- 
mowicz fc Lubow|1994||Hayasaki et al.|2007| MM0 8; [Cuadra 



et al 



2009 



Shi et al.| |2011| 



[Roedig et al. 



2012), for non- 



extreme mass ratios, the binary carves out a cavity in the 
disk, but gas still penetrates the cavity in streams which 
possibly accrete onto the binary components. Here we have 
followed up on the work of MM08 by investigating the nature 
of this inflow across the circumbinary cavity, as a function 
of binary mass ratio q. We have simulated 12 different mass 
ratios in the range 0.01 ^ q ^ 1. This corresponds to the 
expected range of q- values for massive BH binaries produced 
in galaxy-galaxy mergers. 



We find that while the binary 'propellers' are effective 
at maintaining a low-density cavity at the center of the disk, 
they can not efficiently suppress accretion across the cavity. 
For q ~ 1, the accretion rate can even modestly exceed the 
rate for a corresponding single BH. The largest suppression, 
by a factor of ~ 5 compared to a single-BH disk, occurs for 
q ~ 0.05. As long as the circumbinary disk is fueled at a 
near-Eddington rate from large radius, these binaries could 
therefore have quasar-like luminosities. This should facilitate 
finding counterparts to GW events ( Kocsis et al.|2006 \ , and 
should also allow their detection in electromagnetic surveys 
( Haiman et al.|2009 |. 



We have found that the accretion is not only strong, but 
can be strongly variable, with a characteristic (/-dependent 
frequency pattern. While the accretion for q < 0.05 is steady, 
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Figure 12. The time-dependent accretion rates as in Figure[6] except for the single mass ratio q = 0.1 and the four different resolutions 
also shown in Figure [TT] The q = 0.1 binary is at the cusp of the transition from the three-timescale to the single-timescale regime, and 
is particularly sensitive to resolution, as seen especially in the maxima of the accretion spikes. 
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for q > 0.05 there is a strong modulation by the binary, and 
a clear dependence on q of both the variability pattern, and 
the magnitude of the time-averaged accretion rate. For an 
equal mass binary, the accretion rate is modulated at twice 
the orbital frequency and ~ 1/3 the orbital frequency. As 
the mass ratio is lowered, the power in the 1 /2ibm and 3tbm 
variability timescales is reduced, and traded for a third vari- 
ability timescale at fbi n . In the range 0.05 < q < 0.1, only 
the single timescale is present. Interestingly, even when 
the time-averaged accretion rate is suppressed (compared to 
q = 1), we find that the fluctuations increase, and the ac- 
cretion rate during the brief maxima is enhanced. We were 
able to interpret most of the above results based on simpli- 
fied models of the binary torques and kinematics. 

The strong and highly variable accretion, with charac- 
teristic frequencies, should aid in identifying massive BH 
binaries in galactic nuclei. The presence of two frequencies, 
in the ratio 1:2, is an especially robust prediction that is 
independent of disk properties, and could serve as a 'smok- 
ing gun' evidence for the presence of a binary. Our results 
suggest that the ratio of the power at these two frequencies 
could probe the mass ratio q, while other features of the 
periodogram could probe properties of the disk, such as its 
viscosity. 

The variability time-scales are on the order of the or- 
bital period, and we have argued that the most promis- 
ing candidates in a blind electromagnetic search would be 
10 6-7 Mq binaries, with orbital periods of days to weeks. 
The time- variable accretion to the central regions could pro- 
duce corresponding variability in broad-band luminosities, 
allowing such a search in a large time-domain survey, such 
as LSST, without spectroscopy. Additionally, the emission 
lines could exhibit periodic shifts in both amplitude and fre- 
quency; kinematic effects from the binary's orbit could be 
distinguished from those due to the fluctuating accretion 
rate, whenever the latter contains multiples of the binary 
period. 

Parts of the accretion streams generated periodically 
fuel the BHs, but part of the stream material is flung back 
and hits the accretion disk farther out. The shocks produced 
at these impact sites are prominent for q > 0.1, and can 
provide additional observable signatures. In particular, ra- 
diation from these shocks should be temporally correlated 
with the luminosity modulations arising near the secondary 
and/or primary BH, with a delay time on the order of a 
binary orbit. 

GW observatories, such as eLISA, will be able to con- 
strain the mass ratios of in-spiraling MBHB's at the centers 
of galactic nuclei ab-initio, providing a template for the ex- 
pected variability pattern. This should be helpful in identify- 
ing the unique EM counterpart among the many candidates 
(with luminosity variations) in the eLISA error box, as the 
source with a matching period. 

In summary, our results imply that massive BH bina- 
ries can be both bright and exhibit strongly luminosity vari- 
ations, at the factor of several level. This raises the hopes 
that they can be identified in a future, suitably designed 
electromagnetic survey, based on their periodic variability. 
Although encouraging, these conclusions are drawn from 
simplified 2D hydrodynamical models of a real binary disk 
system, and should be confirmed in future work. 
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